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Abstract 

We present the results of using the computer algebra program Cadabra 
to develop Riemann normal coordinate expansions of the metric and 
other geometrical quantities, in particular the geodesic arc-length. All 
of the results are given to sixth-order in the curvature tensor. 

1 Introduction 

In a previous paper [1] we demonstrated, through a series of simple exam- 
ples, how a new computer algebra program, Cadabra ([2], [3J, [I], [3]) could 
be employed to do the kinds of tensor computations often encountered in 
General Relativity. The examples were deliberately chosen to be simple. So 
it is reasonable to wonder how Cadabra might handle much more challenging 
computations. To explore this question we put Cadabra to the task of com- 
puting the Riemann normal coordinate expansions of the metric and other 
geometrical quantities, in particular the geodesic arc length. Such computa- 
tions are known to be, at higher orders, very demanding and prohibitively 
difficult to compute by hand. Here we will show that Cadabra handles these 
computations with ease. 

There is also a direct practical purpose to these calculations. We are de- 
veloping ([BJ, [7j) an approach to numerical relativity that uses a lattice, 
similar to that used in the Regge calculus ([8], [9], [ID]), as way to record 
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the metric (through the table of leg-lengths) and topology (through the con- 
nections between vertices) of the spacetime. Central to that approach is the 
use of Riemann normal coordinates in the computation of the Riemann cur- 
vatures given just the leg-lengths and (some) angles. For this to work we 
need an equation that links the curvatures to the geodesic arc-length (the 
leg-lengths). The result is equation fill. 211) in section (fiTT) . 

The basic idea behind Riemann normal coordinates is to use the geodesies 
through a given point to define the coordinates for nearby points. Let the 
given point be O (this will be the origin of the Riemann normal frame) and 
consider some nearby point P. If P is sufficiently close to O then there exists 
a unique geodesic joining O to P. Let v a be the components of the unit 
tangent vector to this geodesic at O and let s be the geodesic arc length 
measured from O to P. Then the Riemann normal coordinates of P relative 
to O are defined to be x a = sv a . These coordinates are well defined provided 
the geodesies do not cross (which we can always ensure by choosing the 
neighbourhood of O to be sufficiently small). 

One trivial consequence of this definition is that all geodesies through O are 
of the form x a (s) = sa a and that the v a are constant along each geodesic. 
This implies, by direct substitution into the geodesic equation, that T c a b = 
at O which in turn implies that g a b, c — at O. Suppose now that we were 
to expand the metric as a Taylor series in x a about O. In that series there 
would only be the zero, second and higher derivatives of the g a t,. Thus the 
leading terms of the metric could be expressed as a sum of a constant part 
plus a curvature part. If the curvature is weak this can be interpreted as 
an expansion of the metric in powers (and derivatives) of the curvature. 
Likewise one can imagine similar expansions of other geometrical quantities 
(e.g. geodesies, arc length) in terms of a flat space part plus a curvature 
contribution. 

Higher order expansions are extremely tedious to compute by hand. Some 
hardy souls ([H], [12], [IS]) have endured the journey (Miiller et al ([H]) 
venture as far as 11-th order expansions for the metric (i.e. to terms involving 
the 8-th derivative of the curvatures, the notion of order will be defined in 
the following section). However most people settle for just the first two non- 
trivial curvature terms (i.e. R and VR). 

As we will see, the easiest expansion to compute is that for the metric in 
powers of the curvatures and its derivatives. Much more challenging is the 
expression that allows the Riemann normal coordinates to be constructed 
from a given set of coordinates. This later calculation involves the solution 
of a two point boundary value problem - not a job for the faint-hearted. 
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The expansion of the metric in Riemann normal form can be found in many 
articles. For those with a mathematical bent see ([E], [IS], [H]) and in 
particular the elegant exposition by Gray ([IB]. [19]). For applications in 
physics see ([20], [21], [22]). 

2 Conformal coordinates 

Each algorithm given later in this paper yields polynomial approximations 
to particular geometric quantities (e.g. the metric). Higher order approxi- 
mations are obtained by recursive application of the algorithms. 

In this section we will define what we mean when we say that the polynomial 
S e is an expansion of S up to and including terms of order O (e n ). We will 
do so by introducing a conformal transformation of the original metric. 

Consider some neighbourhood of O and let e be a typical length scale for 
O (for example, e might be the length of the largest geodesic that passes 
through O and confined by the neighbourhood). Construct any regular set 
of coordinates x a (i.e. such that the metric components are non-singular) in 
the neighbourhood of O and let the coordinates of O be i". We will use the 
word patch to denote the neighbourhood of O in which these coordinates are 
defined. Now define a new set of coordinates y a by 

x a = xl + ey a 

and thus 

ds 2 = g ab (x)dx a dx b = e 2 g ab (x± + ey)dy a dy b 
Now define the conformal metric ds by ds = ds/e 2 . This leads to 

ds 2 = g a b(x* + ey)dy a dy b = g ab (y, e)dy a dy b 

and 

9ab = 9ab , 9ab,c = ^9ab,c , 9ab,cd = ^9ab,cd at O 

where the partial derivatives on the left are with respect to y and those on 
the right are with respect to x. Since g ab (x+) does not depend on e we have 
the general result that 

9ab,hi2i3-in = O (e n ) at O 
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From this it follows, by simple inspection of the standard equations, that 

f bc,iii2l3— in = O (y6 ) at O 

R a bcd,iii 2 i3—i n = @ ( e " +J ) a ^ ^ 

There are now two ways to look at the patch. We can view it as patch of 
length scale e with a curvature independent of e. Or we can view it as patch 
of fixed size but with a curvature that depends on e (and where the limit 
e — > corresponds to flat space). This later view is useful since in using it 
we can be sure that the series expansions around flat space are convergent 
(for a sufficiently small e). 

We will use these conformal coordinates for the remainder of this paper. As 
there is no longer any reason to distinguish between x a and y a we replace 
y a with x a . The x a will now be treated as generic coordinates (but keep in 
mind that we are working in a conformal gauge). 

Finally, when we say that S e is an expansion of S up to and including terms 
of order O (e n ) we mean that 

\S — S I 
< lim < M 

t^o e n+1 

for some finite positive M i.e. if S were expanded as a Taylor series in e 
around e = then S and S £ would differ by terms proportional to e n+1 . 



3 Riemann Normal Coordinates 



Suppose, as is almost always the case, that our coordinates x a are not in 
Riemann normal form. How might we transform to a local set of Riemann 
normal coordinates? If we were to appeal to the simple definition y a = 
sv a we would soon encounter a hurdle. The quantities v a are rarely known 
explicitly but must instead be computed by solving a two-point boundary 
value problem. This is non-trivial but it can be dealt with in stages. First 
construct a Taylor series expansion for an arbitrary geodesic passing through 
O. This solution to the initial value problem will depend on two integration 
constants, x a and x a , being the respective values of x a (s) and dx a /ds at s = 0. 
Next use a fixed-point iterative scheme to produce successive approximations 
to x a so that the geodesic passes through both points O and P. Then finally 
compute v a as dx a /ds at s — 0. This is exactly the plan which we will follow. 
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It is well known that Riemann normal coordinates can always be constructed 
locally around any non singular point (see [ISj)- Thus we will not concern 
ourselves here with issues such as existence and convergence but rather we 
will focus our attention on developing algorithms for expressing various quan- 
tities in Riemann normal form. 



3.1 The initial value problem 

Our aim here is to obtain a Taylor series, about the point O, for solution of 
the geodesic equation 

q d 2 x a _|_ -pa / \ dx b dx° 
ds 2 c ds ds 

subject to the initial conditions x a (s) = x a and dx a /ds = x a at s = 0. 

We choose s = at O and we write the Taylor series for x a (s) as 



x a {s) = x a \ s=0 + s — 



^ s n d n x a 
„ n + ^ n\ ds n 



s=0 



The second and higher derivatives can be obtained by successive differentia- 
tion of the geodesic equation leading to 



rr a ( q\ — rr a _|_ Q i a _ \ p a . . . . ™H™«2™i3 . . . ~«n (O 1 ^ 

n=2 



where the F a iii 2 i 3 -i n , known as generalised connections, are defined recur- 
sively by 

pa _ -pa ra pp /q n\ 

1 hi2i3—in ~ 1 (hi2i3—in-l,in) lbL p(i2«3---in-i i hin) K ^) 

Note that the use of round brackets (...) denotes total symmetrisation over 
the included indices (see Appendix A for more details). 

A convenient shorthand for equation (13.21) in terms of covariant derivatives 
can be obtained if you ignore (in this context alone) the single upper index. 
This leads to the compact notation 
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3.2 The boundary value problem 



Here we seek to juggle x a so that the geodesic passes through not only O but 
also P. Suppose that the coordinates of P are x a + Ax a where x a are the 
coordinates of O. The solution (13.1 p already passes through O so we have 
only now to force it to pass through P. Let sp be the geodesic distance from 
O to P. Then our challenge is to solve 

A T » _ e „™ a _ \ _Pp» ™H™«2™«3 . . . ™«n 

opu/ 7 ± l;l j 2 j 3 ...j a/ AO/ a. 

^— ' n! 

n=2 

for x a in terms of Ax a and the generalised connections. 

Put y a = spx a (this introduces the Riemann normal coordinates) and re- 
arrange the equation into the form 

y a = Ax a + J2 T a lli2i3 ... ln y i Y 2 y ia ■ ■ ■ V in (3.4) 

n=2 

We now plan to solve this equation for y a by constructing a sequence of 
approximations to y a . In principle we could imagine that y a has been 
found and that we have expanded it as a infinite series in powers of e (i.e. 
as a power series in the curvatures). We will choose y^ to be the Taylor 
polynomial of y a to order e m . That is, y^ is a polynomial in the curvatures 
(and its derivatives) up to and including terms of order 0(e m ). We can 
compute y^ by truncating both side of (13.41) to terms no higher than O (e m ). 
But note that the ^ a i 1 i 2 i z --i n are of order 0(e n ~ 1 ). The upshot is that the 
infinite series may be truncated at n = m while still retaining all terms up 
to and including e m . Thus we have 

(n=m 1 \ 
\^ _F a . . . . 7/S/S/3 ■ • ■ V in I 
n=2 ' / 

Where T e m is a simple truncation operator (it deletes all terms of order 
O (e m+1 ) or higher). This is a marginal improvement on (13. 4p (at least we 
have a finite series) but it is still a non-linear equation for y^. But fortu- 
nately we can do better. Notice, once again, that T a i li2i3 ...i n = O (e n_1 ) and 
this allows us to use lower order estimates for y a in the product terms on the 
right hand side. This leads to 
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Now we see that appears only on the left hand side and thus we can use 
this equation to recursively compute y£ for p — 2, 3, 4, • • -. Here are the first 
few y^. We start with the lowest order approximation, 

y a o = Ax* 

and as there are no e 1 terms in ( 13. we have 

y\ = y a = Ax a 
Now set m = 2 in equation ( 13.51) to obtain 

yl = Ax a + Tl Q^^y?) 

= Ax a + ^T^Ax^Axf 
and once more, with m = 3, with the result 

yl = Ax a + T e 3 (^T\ i2 y^ + ^T^y^y?) 

= Ax a + ir\, 2 Ax ?i Ax' 2 + i (rvr b i2i3 + r° ilia)i ,) Ax Ji Ax J2 Ax i;i 

This process may seem simple but looks can be deceiving - the higher order 
y^ contain a profusion of terms that, when computed by hand, are largely 
unmanageable beyond m ~ 7. We will return to this point later when we 
discuss the use of Cadabra to perform these computations. 

This completes our first objective - to find a way to transform from any 
non-singular set of coordinates to a local set of Riemann normal coordinates. 
The question we now pose is - what form does the metric take in these 
coordinates? This is the subject of the second next section. But first we 
shall take a short moment to introduce some new notation. 



4 Notation 



The calculations we are undertaking are flooded with expression such as 

hhi3—in+l = (iii 2 i3— — + l)r p (i 2 i 3 ...i n T l \ ± i n+1 ) 

in which long stretches of indices like i^iz • • • i n abound. This is tedious to 
write and, in the authors opinion, rather untidy. Here we propose a variation 
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on the notation that is both easy to read and write while not detracting from 
the meaning in the expression. The proposal is that a sequence of indices 
such as ii«2*3 • ■ - in be replaced with a single index of the form i. In this 
notation the previous equation could be written as 

r a 6cd = r a ( fe£ ,d) — (n + i)r a p (cr p bd ) 

where c contains n > indices. In cases where the number of hidden indices 
needs to be made clear we will either say so in words (as we did in the 
previous example) or we will include the number as a subscript, for example 

r a bc n d = r a (6c n ,d) - (n + i)r a p( £n r p fed) 

This version is however, not as clean as the previous version. 
For equations like 

n — r a , ■ ■ ■ ■ A ix A i2 A iz . . . A in 

we will write 

o = r a bc4 A^ 

We chose this notation of single dot in A - because of its suggestive form 
(of multiplication of as many copies of A as required to exhaust the indices 
within d). Here is another common construction 




How might we tidy this up? By including a dot before the derivative index 
d, like this 

These simple changes brings some degree of normalcy to the printed form but 
those gains rapidly pale into insignificance when we ask Cadabra (in section 
(]lip ) to display the results for the sixth order Riemann normal expansions. 
These results contain expressions with symmetries in some of the indices 
which when printed can stretch over many A4-pages. The reason is in part 
that the expressions are inherently long but also because Cadabra does not 
use the round-bracket notation to denote symmetrisation over a set of indices. 
Instead it uses the fully expanded form which, on paper, can lead to an n\ 
explosion in otherwise similar looking terms. We will try to minimise this 
problem as follows. Suppose we have an object A a b c de which we know to 
be symmetric on the indices (cde). We create an arbitrary B a and instruct 
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Cadabra to simplify the expression A a b a i e B c B d B e as much as possible. Then 
we set B a equal to one just prior to printing the expression. But how do we 
convey to the reader that this operation has been performed? The expression 
that we are trying to print will always be the right hand side of some equation 
such as D a b c de = ^4ab(cde) and thus D a b c de is also symmetric on (cde). So 
when we print this equation we will include the round brackets on the left 
hand side. Thus we would display the results as D ah ^ cde ) = A abcde and we 
would understand that the right hand side should be symmetrised over (cde). 
Clearly this notational device should only ever be used at the end of the 
calculations. 

Our final notational device concerns cases where we want to exclude an index 
from symmetrisation. The normal practise is to exclude an index by enclosing 
it in a pair of vertical lines. In our variation we will place a dot above 
the excluded index. Thus, where other authors might write (ab\c\d\e\fg) to 
signify symmetrisation over only a, b, d, f and g, we will write (abcdefg). 

Though these variations might appear to make only marginal improvements 
in the printed equations they have made it considerably easier for the author 
in creating the LaTeX code for this document. 

5 The metric in Riemann normal form 

In the preceding section we chose to distinguish between generic and Riemann 
normal coordinates by using the symbols x a and y a respectively. We will now, 
for notational convenience and to accord with convention, revert to using 
x a for the Riemann normal coordinates while stripping y a of any special 
meaning. 

Our aim in this section is to express the metric in Riemann normal form. 
This will take the form of an infinite series in powers of the curvature tensor 
and its derivatives. We start by writing out the Taylor series for the metric 
around x a = 



where c contains n indices and g a b are constants (e.g. g a b = diag(l, 1, 1, • • •)). 

Our present task is to express the partial derivatives of the metric in terms 
of the Riemann tensor. From the standard definition of a metric compatible 
connection we have 




n=l 
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and since g a b,cd is totally symmetric in cd we also have 

9ab,cd = (gae^ e b(c) ,d) + (deb^^c) ,d) 

Two points should be noted, first, the connection appears only in the form 
r a b(c,d), second, the left hand side contains derivatives one order higher than 
in the corresponding terms on the right hand side. The upshot is that we can 
use this equation to recursively compute all of the metric derivatives solely 
in terms of the T a b ^ C)d ) and the constants g^. In this way we could express 
the above Taylor series for the metric in terms of the connection. But we 
can go one stage further - the derivatives of the connection must surely tie 
in with the curvatures. Thus we are led to review the standard definition for 
the curvature, which after a series of derivatives, can be written in the form 

j~ta _ pa pa i / pa pi \ / pa pi \ 

K {bcdS ~ d( - bc 'S) ~ ( 6c >£) d + V i( - c bd > '£) ~~ V id ( - bc ) '£) 

Can we use this to eliminate the connection and its derivatives from the 
metric? Yes, but only after we specialise to the Riemann normal coordinates. 

Recall that, in Riemann normal coordinates, all geodesies through O are of 
the form 

x a ( s ) = sv a 

which upon substitution into the geodesic equations leads to 

= T a {b£) at O 

It follows, by recursion on equation (13.31) . that 

= T\ bcS at O (5.1) 

We also know that T a b c , e d is separately symmetric in its first pair of indices 
and in the remaining (n + 1) lower indices (assuming e contains n indices). 
Thus using equation (1 A. 1 1) we have 

= (n + 3)r a ( bCi g d ) = 2T a d ^ C e) + (n + l)r a ( feC) e)d 

We can use this to eliminate the r a ( &c e ) rf term in the previous equation for 
the curvature. The result, after a minor shuffling of terms is 

(n + 3)r = (n + 1) (R\ bcdS - {T a i{c T\ d ) >eJ + {T a td T\ bc ) s 

(the reason for rearranging the terms will become clear in a moment). Note 
also that the last term in the previous equation can be eliminated by equation 
(15.11) and the product rule. 
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In summary, the equations of interest are 

9ab(x) = 9ab + ^2 — 9ab,c X'~ (5.2) 
n=l U ' 

9ab,cd (9ae^ e b(c) ,d) + {g e b^ e a(c) ,d) (5-3) 

(n + 3)r>, c ,) = (n + 1) (^ Mg) - (rVJ ,e)) (5.4) 

We use these equations as follows. First, equation f)5.4p is used to recursively 
compute the T a b( c ,ed) i n terms of the Riemann tensor and its partial deriva- 
tives (this was the reason behind the shuffling of terms noted above). Note 
that e in equation (15.41) contains n hidden indices. The T a ^ c ^) are then sub- 
stituted into (15.3p which in turn is used to recursively express all of the g a b, c 
in terms of the Riemann tensor and its partial derivatives. When the dust 
settles we will have a finite series expansion for the metric in terms of the 
Riemann tensor and its partial derivatives. The result, accurate to O (e ), is 

gab (x) = g a b- I x c x d R acbd - - x c x d x e d c R adbe + O (e 4 ) 

6 

We conclude this section by introducing one final variation to the algorithm 
just given - an option to re-express the metric in terms of the covariant rather 
than partial derivatives of the curvatures. 

It is not hard to see that after a series of covariant derivatives one would 
obtain an equation of the form, in any coordinate frame, 

pa pa I /Qa 

(bcd;e) ~ ^ (bcd,e) ' V (bcde) 

where Q a {j bcde ^ j is a function of the T a bc , the R a b c d and their partial derivatives. 
If this is going to sit nicely with our algorithm given above then we will need 
to show, in the Riemann normal frame, that this equation only contains 
connection terms of the form T p q ^ r>s y Fortunately this is rather easy to do. 

Each term of the form T p qrtS in Q arose during one round of covariant dif- 
ferentiation. Thus at least one of the indices q, r and all of the indices in 
s must be drawn from the index list e. If both q and r are contained in e 
then the term is of the form r p (q rjS ) and thus will vanish when we specialise 
to the Riemann normal frame. This completes the proof. If we re-arrange 
the above equation into the following form 

R a {bcd,e) = ^(bcdie) ~ (bcde) (5-5) 
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we can use it to recursively compute all of the partial derivatives of the cur- 
vatures in terms of their covariant derivatives. The Q a ^ bcde - ) will contain lower 
order derivatives of the curvatures and partial derivatives of the connection 
all of which can be eliminated (in favour of covariant derivatives) using previ- 
ously computed results. For the first two derivatives we find that the partial 
and covariant derivatives are equal but differences do appear in higher order 
derivatives. More details will be given in a later section. 

These calculations, as simple as they may appear, are exceedingly tedious to 
do except for the first few terms. The recursive nature of the calculations 
requires frequent substitution of one result into another which causes an 
explosion in the number of terms that must be handled. Not only is this 
tedious but its is also extremely prone to human error. Calculations of this 
kind are clearly best left to a computer. We shall return to this point later 
on. 



6 The inverse metric in Riemann normal form 

Most of the hard work is now behind us and we can now develop algorithms 
for Riemann normal expansions for other interesting quantities, in this in- 
stance the inverse metric g ab (x). In the previous section we used = g ab - c as 
our starting point. On this occasion we start with = g ab - c - Then, following 
a path similar to that used in the previous section, we arrive at the following 
equations 

oo 1 

9 ab (x) = g ab + V -g a \ c x- c - (6.1) 

n=l 

g a \ cd _ = - {g ae r b e{c ) s - {g eb r a e{c ) s (6.2) 

(n + 3)r d(6) ce) = (n + 1) (R\ bcdS - (r\ c T\ d ) ,,)) flBZJ) 

These equations can be used to construct the series expansion for g ab (x), 
which to O (e 4 ) is 

g * b ( x ) =g ab + \ x c x d R a c b d + l - x c x d x e d c R a d \ + (e 4 ) 
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7 Generalised connections 



In section l3~T1 we saw that the generalised connections T a bcd arose from succes- 
sive differentiation of the geodesic equation and that they can be computed 
recursively using 

rv = v \bc,d) ~(n + i)T a p(£ v p bd) dS2D 

where the index list c contains n > indices, or directly using 

r a fe c d = r a ( fec; d) (I3.3P 
Here are the first three generalised connections 

F(bc) ( x ) — r a be 
^(bcd)( x ) = dbT a cd — 2 r a be r e cd 

r? bcde) (x) = -T f hc d f T a de -4T f bc d d T a ef + d bc T a de 

+ 2 T a fg Tf bc T 9 de + 4 T a bfTf cg T° de -2T a bf d c T f de 

and when we specialise to Riemann normal coordinates we obtain 

rffc) ( x ) = l xdRa bde + ^ x d x e (2 V b R a dec + 4 V d R a bee + V" R dbec ) + O (e 4 ) 

T a {bcd) (x) = ±x e V b R a ced + 0{e 4 ) 
T a {bcde) (x) = O {e*) 

8 Geodesies 

When first we spoke of Riemann normal coordinates we restricted our at- 
tention to the geodesies that passed through the point O. Now we wish to 
be somewhat less restrictive. We would like to know how to construct any 
geodesic in the neighbourhood of O. Here we will once again be building 
solutions of the geodesic equations and as before we will consider two sepa- 
rate cases, first, the geodesic initial value problem and second, the geodesic 
boundary value problems. 

Most of the machinery that we need to tackle these questions has already 
been developed. Here we apply the formalism developed in sections 13.11 and 
13.21 to the metric in Riemann normal form as obtained in section [51 
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8.1 Geodesic initial value problem 



Consider a point P distinct from O. At P we can assume that the generalised 
connections do not vanish (which is generally true, the exception being flat 
space) . Thus the coordinates x a in the neighbourhood of P do not constitute 
a Riemann normal frame relative to P. But as P lies in the patch for O we 
know that the metric is non-singular at P and thus we should be able to 
construct a new set of Riemann normal coordinates, y a , with P as the origin. 

We have seen this problem once before, in section 13.21 Using equation (18. ip 
and the generalised connections from section [7J we find 



x a (s) = x a + sx a - — s 2 x b x c (8 x d R a bdc + 2 x d x e \7 b R a dec + 4 x d x e V d R a bec 



8.2 Geodesic boundary value problem 

Consider now the case where we have three distinct points O, P and Q. In 
this section we seek to compute the geodesic that passes through P and Q. 
We use the equation (I3.5P of the generalised connections from section [7J to 
obtain 



1 



+ x d x e V a R dbec ) 



l -s i x b x c x d x e V b R a ced + 0(e A ) 



x a (s) =x a + sxl + s 2 x a 2 + s 3 4 + O (e 4 ) 




Ax a + - x b Ax c Ax d R a cbd + — x b x c Ax d Ax e V d R a ba 

O J. Zi 

+ \ x b x c Ax d Ax e V b R a dce + ^ x b x c Ax d Ax e V a R bdce 
o 24 

+ — x b Ax c Ax d Ax e V c R a dbe 
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V b R a dce - 7^7 x b x c Ax d Ax e V a R 



bdce 



x 



b Ax c Ax d Ax e V c R 



•a 



dbe 
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9 Geodesic arc-length 



Since we now have explicit expressions for the metric and the geodesic that 
joins the points P and Q we can compute the length of that geodesic by way 
of the integral 



f Q ( , ,dx a dx b \ 1 , 
L PQ = J p (**(*) — —j ds 



We have, to this point, taken s to be the proper distance along the geodesic. 
However, after careful inspection of the geodesic path ( 13.1 ft we see that any 
uniform scaling of s is allowed. Thus we can re-scale s so that s = at P and 
s — 1 at Q (of course, the parameter s no longer measures proper distance). 
Furthermore we know that the integrand is constant along the geodesic and 
can thus be evaluated at any point which we shall chose to be P. Thus the 
integral is trivial and we have 

L PQ = 9a b (x) — — 

Using our previous results we obtain 

Ll Q = S .Ax • Ax ' - i A'A, 'Ax ^ - ± ^Ax 'Ax ^Ax • VA* 
-Ix-.xVAx^Ax'V^ H „ + 0( e «) 



10 Cadabra 

No matter how determined one might be, the recursive nature of the preced- 
ing equations will wreak havoc with one's sanity should one dare to venture 
beyond the first few terms. For the good of all concerned it is far better 
to leave such calculations to computer programs such as Cadabra. In this 
section we discuss some particular issues we encountered when writing our 
Cadabra programs. 

The examples given in the previous paper [1] actually arose from our earlier 
investigations of Cadabra as a tool to compute Riemann normal expansions. 
Thus it is no surprise that the techniques given in that paper are well suited 
to (most of) our current needs but with two exceptions. First, we need to 
develop Cadabra code for the truncation operator T e m introduced in section 
( 13.21) . Second, we need to extend the ideas given in pQ to allow Cadabra to 
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compute the symmetrised derivatives of arbitrary tensors, such as R a ^ C( i- e )- 
We shall deal with these issues in the following sections after which we will 
present the final Cadabra generated expressions (to O (e 6 ) and in all their 
gory detail). 



10.1 Truncation of polynomials 



In section H3. 21) we introduced a truncation operator T e m . Here we discuss 
how we constructed that operator in our Cadabra code. It happens to be 
an extremely useful piece of code and is used extensively throughout our 
Cadabra code. 

Suppose you are asked to extract the leading terms from an expression such 

as 

p*(x) = c a + c a b x b + c a bc x b x c + c a bcd x b x c x d 

A standard approach would be to compute the derivatives of P a (x) at x — 
0. This approach would be rather simple to code in Cadabra. However a 
minor issue does pop up. Unless otherwise told, Cadabra will assume that 
all objects have non-zero derivatives. Thus if Cadabra were instructed to 
compute d/ds of the above expression it would dutifully do so but it would 
treat the x's and the coefficients c b as possibly depending on s. Cadabra can 
be coaxed to restrict the derivative operators to act only on specific objects 
by using the : : Depends property and the Ounwrap algorithm. 

A better solution, one that does not involve derivatives, is to use Cadabra's 
: : Weight property and the @keep_weight algorithm. The idea is to assign 
weights to nominated objects (through the : : Weight property) and then 
extract terms matching a chosen weight (using the @keep_weight algorithm). 
Here is a small piece of Cadabra code that does the job. 

x~{a}: :Weight(label=xterms,value=l) . 



poly:= c~{aj- 

+ c~{a}_{b} x~{b} 

+ c~{a}_{b c> x~{b} x"{c> 

+ c~{a}_{b c d> x~{b} x"{c} x~{d>; 



term00:=@(poly) : @keep_weight ! (termOO) {xterms}{0} 
termOl : =@(poly) : @keep_weight ! (termOl) {xtermsMl} 
term02 : =@ (poly) : @keep_weight ! (term02) {xterms}{2} 
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The first line identifies the x a terms as our target (which we name as xterms 
so that they can be distinguished from other targets decalared by other in- 
stances of : : Weight). The next three lines then extract the th , 1 st and 
2 nd terms in the expression poly. The result would be exactly as if we had 
written 

termOO : =c"{a} ; 
term01:=c"{a}_{b> x"{b>; 
term02:=c~{a}_{b c> x~{b} x"{c>; 

The truncated polynomial (in this quadratic) could then be computed 

simply as @(term00)+@(term01)+@(term02). 



10.2 Symmetrised covariant derivatives 

Let v a be a tensor field and suppose we wish to compute v a -b at the origin 
of the Riemann normal frame, O. Construct any geodesic through P and an 
auxiliary field A a throughout the patch containing O. Let the geodesic be 
parametrised by the proper distance s and described by x a = x a (s) with unit 
tangent vector D a = dx a /ds. We choose A a so that it is auto-parallel along 
the geodesic. Thus we have 



dv a 
ds 



dD a 

= V D D a = —— + T a bc D b D c 
ds 

A Aa 

= \7 D A a = — + T a bc A b D c 
ds 

from which it follows that 

Va . b A a D b - = at O 

'- ds n 

where b contains n indices. Thus any higher order covariant derivative can 
be obtained simply by expanding the right hand side, one derivative at a 
time, while using the parallel transport conditions listed above to eliminate 
derivatives in A a and D a . The Cadabra code for this is much the same. Each 
successive covariant derivative is obtained by applying d/ds to the previous 
result then using substitutions to eliminate the newly introduced derivatives 
of A a and D a . Comparing coefficients of A a D- across the equals sign will 
then reveal an expression for v a (-b), the fully symetrised covariant derivative 
of v a at O. 
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Note that = T a ( bc ^ at O in Riemann normal coordinates and thus = 
d n D a /ds n . This can be used to considerably simplify the above computa- 
tions. 

This idea can be easily extended to other cases. For example, suppose we 
require v ab ( c - de ) at O. In this case we would introduce two auxiliary fields, 
A a and B a , each constrained to be auto-parallel along the chosen geodesic. 
Then following the above procedure we obtain 

v ab{c ., de) A«B b D'D*De = gq*^gg) at o 



10.3 Symmetrised derivatives of the Riemann tensor 

By inspection of equation (15.41) we see that the only derivatives of the Rie- 
mann tensor that enter our calculations are of the form R a (bcd,e) and R a ( bcd - e )- 
According to the prescription given above we can compute the covariant 
derivatives in terms of the partial derivatives by the following procedure. 



= V D D a 



dD a 
ds 



= S7 D B\ = ^ + T a dc B d b D c - T d bc B a d D c 
ds 

d n 

R\ bcd .^B c a D b D d D e - = — (R\ cd B c a D b D d ) 

This is not exactly what we want - it yields the covariant derivatives in terms 
of the partial derivatives. We need instead, the partial derivatives expressed 
in terms of the covariant derivatives. In section [5] we provided one solution 
to this problem. There we argued that the equations could be re-written in 
the form 

R a (bcd,e) = R a (bcd;e) ~ Q" (bcde) 

where Q a ^ bcde ^ contains all the lower order symmetrised partial derivatives of 
R a bcd . This algorithm certainly works but it is computationally expensive. If 
for the moment we label the equation for the n-th derivative as E n then our 
algorithm entails a whole hierarchy of substitutions of Ej into Ek for j < k. 
For example, to compute E A we need to substitute Ei into E 2 , then E\ and 
E 2 into E 3 , and finally, E\, E 2 and E 3 into E4. Our code took about 4 seconds 
to compute the first four derivatives but with one extra derivative this time 
grew to 11 minutes. We did not bother to compute the sixth derivative (in 
fact for the results given here we only need the first three derivatives). 
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It seems reasonable to ask : is there a better way? Indeed there is and the 
changes required are very simple. 

The auxiliary field B a b can be freely chosen so there is nothing stopping us 
from setting B a b to be constant throughout the neighbourhood of O. Thus 
every partial derivative of B a b is zero at O. What use is this to us? Consider 
the equation 

d n 

(R\ cd B c a D b D d ) , e _D e - = — (R\ cd B c a D b D d ) 

This is clearly true for any B a b and D a = dx a /ds (and in any coordi- 
nate frame). In our Riemann normal frame we have = d n D a /ds n and 
we have chosen = B a ce . Thus the right hand side can be reduced to 
R a bcd,eB c a D b D d D-. This leads to the following modified scheme (after swap- 
ping the left and right hand sides) 

dD a 

= D a :b D b 



ds 

D D b = D b ;c -Ly = l^tfc-D b 1 ^ — 1 be -D d-> 

R a bcdte _B c a D b D d D e - = (R a bcd B c a ) ]e _D b D d D^ 



Vr>a 75a n c F a n c R a 7~) c 



This algorithm computes all the partial derivatives directly, without requiring 
any substitutions from previous results. It took less than 10 seconds to 
compute the first five derivatives which is a dramatic improvement over our 
previous algorithm (11 minutes). 

In our Cadabra code we compute all of the symmetrised partial derivatives 
(using the above algorithm) before we compute the metric expansion (equa- 
tions (15.21 15.31 15.41) ). In this way we obtain a series expansions for the metric 
in terms of the covariant derivatives of the curvatures. 



11 Expansions to sixth order 

All of our O (e 6 ) Cadabra programs were not overly demanding on compu- 
tational resources, taking less than 2 minutes to run (on a MacOSX with an 
Intel cpu) and requiring less than 13 Mbyte of memory. 

The Cadabra codes and several support scripts are available from the author's 



web site http://users.monash.edu.au/~leo 



The metric 
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180# ai , 0) = 180 g ab - 60 x c x d R acbd - 30 x c x d x e V c R adbe + 8 x c x d x e x f R acgd R begf 
-9x c x d x e x f \/ cd Raebf + 4 x c x d x e x f x 9 R achd V e R bfhg 
+ Ax c x d x e x f x 9 R bchd V e R afhg -2x c x d x e x f x 9 V cde R afbg + O (e 6 ) 

(11.1) 

The inverse metric 

180^ (x) = 180 g ab + 60 x c x d R a c b d + 30 x c x d x e V c R a d b e 

+ 12 x c x d x e x f R a cdg R b efg + 9 x c x d x e x f V cd R a e b s 

+ 6 x c x d x e x f x 9 R a c ^V e i? 6 fgh + 6 x c xVx Vi? 6 crf/i V e i? a fgh 

+ 2 x c xVxV V C(feJ R a / b g + O (e 6 ) 

(11.2) 

Generalised connections 

180r^ c) (x) = 120 x d R a bdc + 15 xV (2 V b R a dec + 4 V d R a bec + V a R dbec ) 

+ X d X e X-F (32 i? a degRfhcg — 16 i? Q bdg R ec fg ~ 8 -R dbg Reef g 

+ 18 V db R a e fc + 18 V de R a b f c — 8 R a gdb R ec fg + 9V a d R ebfc ) 
+ xVx V (16 R dbch V e R a fgh + 6R a deh V b R fcgh 

+ 16 i? Q de hS jRgbch — 8 R a bd hV e Rf C gh ~ 4 -R° d bhS eRfcgh 

— 4 R dbe hV c R a fgh — 8 R dbe hS fR a cgh ~ 4 R dbe hVfR a g C h 

+ 6 V deb R a f gc + 4 V de fR a bgc — 5 R a de h^h,Rjbgc 

— 4 i? a hdb^ eRfcgh — 4 R d behS a Rfcgh ~ 4 R dbe hS /-R" hgc 

+ 3V a de R fbgc ) + 0(e 6 ) 
(11.3) 

180r^ d) (x) = 90x e V feJ R a ced 

+ 3 X 6 X^ (8 i? a ebgRfcdg + 32 i? a begRfcdg ~ 8 R a bcgRedfg 

+ 18 V e fe-R a c/d + 6 V bc R a e fd + 24 i? a gebRfcdg + 3V" bRecfd) 

+ 10 x e x V (2 R ebch V d R a fgh + 2 R ebch V h R a fgd 

+ 4 R ebc hV fR a dg h + 4 R ebc hVfR a h g d + 2 R e bchS a Rfd g h 

+ 2 R a be hS cRfdgh + 4 i?" behS fRgcdh ~ -R" behS hRfcgd 
+ 2 R a heb^cRfdgh + 4 i?" heb^fRgcdh ~ -R° ZieftV hRfcgd) + (e 6 ) 
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(11.4) 

18r" fec(fe) (x) = 8 X S R a bcgRfdeg 

+ X f X 9 (2 R a bch V ' dRfegh + 4 R a bch^fRgdeh ~ bchS? hRfdge 

+ 2 RfbchS dR a geh + 10 RfbchSdR a hge + 4 RfbchS g-R" deh 
+ 8 Rfbch^hR a dge + 2 RfbchS° ' Rgdeh + 12 -R/ftc/i Vrf-R Q epft 
+ 6 -R" bfhS cRgdeh + 6 -R" hfb^ cRgdeh) + C (f 6 ) 

(11.5) 

3r^ e/) (x) = (2 iT bchVdRgefh + 3 RgbchV d R a efh ) + O (e 6 ) (11.6) 

r^ de/!? )^) = 0(e 6 ) (ii.7) 

Partial derivatives of the Riemann curvature tensor. 

R U {bcv,a) = V a R U bcv (H-8) 
R U {cdv,ab) = ^ abR U cdv (H-9) 
2R U (dev,abc) = 2 V sfci?" d e „ — R va bfV ' C R U def + -R" abfV c R vde f (11.10) 
5-R (efv,abcd) = ^^abcdR efv 7 -Rua&g V c <2-R efg + 7 -R a bg^ cdRyefg 

(11.11) 

Riemann normal coordinates 

y a = Ax a + y a bc Ax b Ax c + y a bcd Ax b Ax c Ax d + y a bcde Ax b Ax c Ax d Ax e (11.12) 
+ y a bcdeJ Ax b Ax c Ax d Ax e Ax s + O (e 6 ) 

2C = T% C (11.13) 

Qy a bcd = T a be T\ d + d b T a cd (11.14) 

2% 6 a Cffe = 2 T a bf d c T f de + T a fg T f bc T 9 de + T f hc dfT a de + d bc T a de (11.15) 
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360^/ = -4 r a bg T° ch r h dl r e/ + 2 r a b9 r» ch d d r h e/ + 3 v a bg r» hi r h cd r ef 

+ 6 T a bg T h cd d h T* ef -QT a bg T h cd d e T° fh + 9T a bg d cd T° ef 

+ 13 T a gh T 9 bc d d T ef 

+ r» bc v h dg d h v a ef - 4 r» fcc r h ^d e r a /fc + 7 «9 3 r a ^r* e/ 
+ 2 <9,r a Cff c> d r 9 e/ + 3 bc r h de d g r a fh + 3r° bc r h de d f r a gh 
+ 6 r 9 bc d dg r a ef - 3 r 9 6c <9 de r a /g + 3 d bcd v a ef 

(11.16) 

360^ = 120 x d R a bdc + 15 x d x e (2 V fe i? a dec + 4 V d i? a bec + V a i^ ec ) 
+ x d x e x^ (32 _R a dge R g bfc — 16 i? a bgd R ge f c — 8 R a dgb R ge fc 

+ 18 V dbR a efc + 18 V deR a bfc + SR a gdbRgefc + 9V" d R e bfc) 

+ x d x e x f x 9 (16 R hbdc V e R a fhg + 6R a dhe V b R hfgc + 16 R a dhe V f R hbgc 

— 8 R a bhdS eRhfgc ~ 4 R dhb^ eRhfgc ~ 4 Rhdeb^ cR° fhg 
— 8 RhdebVfR a chg ~ 4 Rhdeb^fR " g hc + 6 V de bR a fgc + 4 V de fR a bgc 
+ 5 R a dhe^hRfbgc + 4 _R Q hdb^ eRhfgc — 4 Rhdeb^"" Rhfgc 
+ 4 Rhdeb^ fR a hgc + 3 V" de Rfbgc) 

(11.17) 

1080y^ cd = 90 X £ VbR a ced + 3 X e X'^ (8 i? a egbRgcfd + 32 i? a bgeRgcfd — 8 i? a bgcRgefd 
+ 18 V eb R a c f d + 6 V bc R a efd + 56 i? a gebRgcfd + 3 V a bRecfd) 
+ 10 X 6 X V (2 R hb ecV d R a fhg + 4 RhbecV h R a fgd + 4 i4fe eC V /iT 

+ 8 Rh bec ^7fR a hgd ~ Rhbe<N° ' Rhf gd + 2 _R Q bhe \ ' c Rhfgd 
+ 4 -R° bhe^ fRhcgd + bhe^ hRfcgd + 4 -R° heb^cRhfgd 
+ 8R a heb^ fRhcgd + 2 _R a heb^ hRfcgd) 

(11.18) 

432< cd£ = 8 s'iF fcgcJ R 9d/e + xV (2 bhc V d R h fge + 4 i? a bh cVfR hd ge 

+ -R bhc^hRfdge + 2 Rhbfc^dR ghe ~ 10 Rhbfc^7 d R hge 
+ 4 Rhbfc^ gR a dhe + 28 -R/ifo/ c Vft-R" dc/e + 2 Rhbfc^ a Rhdge 
+ 12 RhbfcVdR a ehg + Q R a bhfV c Rhdge + 18 -R" ft/fe V c Rhdge) 

(11.19) 

360C rfe/ = (2 ^ bhNdRhegf + 3 R hb gcV d R a efc/) (11.20) 
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Geodesic arc-length 



L 2 PQ = f ab Ax a Ax b + f abc Ax a Ax b Ax c + f abcd Ax a Ax b Ax c Ax d (11.21) 
+ f abcde Ax a Ax b Ax c Ax d Ax e + O (e 6 ) 

180/ a6 = 180 g ab - mx c x d R cadb - 30 x c x d x e V c R daeb 

+ X C X d X £ X f (8 RgcdaRgefb ~ 9 V cd R ea f b ) 

+ 2x c x d x e x f x 9 {AR hcda V e R hfgb - V cde R fagb ) 

(11.22) 

180/afec = -15 X d X e V a Rdbec + X*X e X f (8 RgdeaRgbfc ~ 9 V ' daRebfc) 

+ X d X E X f X 9 (4 Rhadb^eRhfgc + 4 RhdaN bRhfgc + 4 RfideaS fRhbgc 

— 3 V 'deaRfbgc) 

(11.23) 

540/ afecd = -3 (44 RgaebRgcfd + 3 V ab R ecfd ) - 5 X e X J X 9 (8 Rhaeb^cRhfgd 

+ 9 Rhaeb^ \Rfcgd + 20 Rhaeb^ fRhcgd ~ 6 Rhe f bRhcgd) 

(11.24) 

54/ abcde = X f X 9 R hafb V c R hd ge (11.25) 

12 Discussion 

The value of any new computational tool comes not just in being able to do 
routine computations, computations that we could do by hand, but rather 
in giving us the option to perform computations we would not otherwise 
undertake. New tools should open new opportunities for research. Cadabra 
seems to be such a tool. 



13 Source 

A .tar.gz archive of the Cadabra files used in preparing this paper can be 



found at this URL http : //users .monash. edu. au/~leo/research/papers/files/lcb09-03 .html 
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Appendix A. Symmetrisation of tensors 

The totally symmetric part of a tensor Ai x i 2 i 3 ...i n is commonly defined by 

A(i 1 i 2 i 3 ---i„) = (^4iii2«3"-«n "I" Ai 1 i 2 i a ...i n + Ai t i 2 i 3 ... i n + • • •) 

where the sum on the right hand side includes every permutation of the 
indices of ■ ■ • i n - If the tensor A ili2 i 3 ... in happens to be symmetric in 
every pair of indices then we observe 

A,. . . A. . . . 

From the above definition it is very easy to establish the following theorems 

A{iti % i z — (313233— jm)— in) = A(i 1 i 2 i 3 ---j 1 j 2 j 3 ---j m ---i n ) 

A(i 1 i 2 i 3 ---i n Bj 1 j 2 j 3 ...j m ) = A((i 1 i 2 i 3 ...i n }B(j 1 j 2 j 3 ...j m ^ 

n A(i 1 i 2 i 3 i 4 ...i n ) — ^4ji(i 2 i3i4---in) ~l~ Ai 2 (i 1 i 3 i 4 ...i n ) + Ai 3 (i 1 i 2 i 4 ...i n j + • • • 

~l~ Ai n (i 1 i 2 i 3 ...i n _ 1 ) 

^-^-(11121314— in) A(i 2 i 3 i i ...i n y ll + -*4(iii3i4..-i n )i2 "I" -"-(ii i 2 u---i n )i 3 + ' ' ' 

+ 7 4(jlj2«3---*u-l)«n 

Suppose now that we have A ili2 i 3 ... in = A {hi2i3 ... in) , that is, A ili2 i 3 ... in is totally 
symmetric. Then for any Bj we have 

(n + l)A(i 1 i 2 i 3 ...i n Bj^ — Aji 2 i 3 ...i n Bi x + Ai x ji 3 ..., ln Bi 2 + A^ x ^ 2 j...^ n B^ 3 

+ ■ ■ • + Ai x i 2 i 3 ...i n _^Bi n 

and 

(n + 1)74.(2^223. ..j n j) — Aji 2 i 3 ...i n ^ x + Ajjjjg...^^ + y4j 1 j 2 j...i ni 23 

+ • ■ ■ + Ai 1 i 2 i 3 ...i n _ 1 j ! i n + Ai 1 i 2 i 3 „,i n ^ 
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All of the above are very easy to prove but one result which requires just a 
little more thought is the following. 

Suppose Ai l i 2 j. i j 4 j 5 ...j n is symmetric in the pair i\%2 and symmetric in all the 
indices J3J4J5 • • • j n - That is, it is symmetric under the interchange of any 
pair of z's and any pair of j's but it is not necessarily symmetric when any i 
is swapped with any j. What can we say about A(i li;2 j 3 j 4 j s ...j n \l Here is the 
result 

n A(i l i 2 i 3 —i n ) = 2^4i n (iii 2 i3— in-i) ~\~ { n ~ tyAfai^. ■■i n _ 1 )i„ (A-l) 

The proof is very easy. Begin by writing out n\ ^4(j 1 i 2 j 3 - -i n ) in full. Then 
partition the terms into two disjoint sets, one set in which i n appears in 
one of the first two index slots, the other set in which i n appears in any of 
the remaining n — 2 slots. The terms in the first set are exactly those that 
define A^^...^) while those in the second set define A^ li2 i 3 ... in _^ in . The 
above equation follows by simply counting the number of terms in each set 
(2(n — 1)! and (n — 2)(n — 1)! respectively) and the simple observation that 
n\ A( ili2 i 3 ... in ) equals the sum of the terms from both sets. 

Finally we note that if Q = A ili2i:j ... in x ll x l2 x 13 ■ ■ -x %n then we have 

Q ,hhi3---in ~ ^' ^(iii2i3 f "i») 

Cl — A, x T* 1 T* 2 T* 3 . . • T in 
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